#' =============================================================================
#' FILE: 03_cueing_analysis.R
#' DESCRIPTION:
#'   Analyzes the effects of ingroup/outgroup cues on policy support using
#'   regression and hierarchical models. Produces main tables and figures for
#'   the cueing experiment.
#'
#' PACKAGES REQUIRED: pacman, tidyverse, miceadds, lme4, stargazer, ggpubr
#'
#' OUTPUTS:
#'   - 03_figures/fg_02.jpg
#'   - 03_figures/fg_04.jpg
#'   - 03_figures/fg_A01.jpg
#'   - 03_figures/fg_05.jpg
#'   - 04_outputs/table_a04.doc
#'   - 04_outputs/table_a04_complete.doc
#'   - 04_outputs/table_a05.doc
#'   - 04_outputs/table_a05_complete.doc
#' =============================================================================

# Packages and dataset ---------------------------------------------------------

# Required packages
if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, 
               miceadds,
               lme4,
               stargazer, 
               ggpubr)

# Data
df <- readRDS("04_outputs/clean_dataset.rds")
df_e <- readRDS("04_outputs/clean_dataset_reshaped.rds")

# Variables used in the models -------------------------------------------------

# Ingroup vs Outgroup cue
df_e$inout <- ifelse((df_e$cha_id == "Fujimorismo" &
                        df_e$cond == 'fuji') |
                       (df_e$cha_id == "Anti-Fujimorismo" &
                          df_e$cond == 'anti'), "ingroup", 
                     ifelse((df_e$cha_id == "Fujimorismo" &
                               df_e$cond == 'anti') |
                              (df_e$cha_id == "Anti-Fujimorismo" &
                                 df_e$cond == 'fuji'), "outgroup", 'other'))

df_e$inout <- factor(df_e$inout, levels = c('other', 
                                            "ingroup", 
                                            "outgroup"))

# In and outgroup cues together
df_e$inout2 <- ifelse(df_e$inout == 'other', 
                      'other', 
                      'inout')
df_e$inout2 <- factor(df_e$inout2, levels = c('other', 
                                              'inout'))

# Condition
df_e$cond <- factor(df_e$cond, levels = c("none", 
                                          "fuji", 
                                          'anti'))
# Movement-based Identities
df_e$cha_id <- factor(df_e$cha_id, levels = c("none", 
                                              "Fujimorismo",
                                              "Anti-Fujimorismo"))

df_e$exp_type <- ifelse(df_e$exp == 1 |
                          df_e$exp == 4 |
                          df_e$exp == 6, "fuji", "anti-fuji")

# Strength
df_e$fuji_st <- ifelse(df_e$cha_id == "Fujimorismo", df_e$cha_st, 0)
df_e$anti_st <- ifelse(df_e$cha_id == "Anti-Fujimorismo", df_e$cha_st, 0)

# Figure 02: Fujimorismo and anti-Fujimorismo in Peru --------------------------

fuji_plot <- df %>% 
  select(cha_id, 
         cha_st) %>% 
  pivot_longer(-cha_id) %>% 
  filter(cha_id != 'none') %>% 
  mutate(cha_id = factor(cha_id, levels = c("Fujimorismo", 
                                            "Anti-Fujimorismo"))) %>% 
  ggplot(aes(x = value, color = cha_id)) +
  geom_density() +
  labs(x = 'Identity Strength', y = "Density") +
  scale_color_manual(values = c("Black", "Gray")) +
  theme(legend.position = "bottom", 
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5),
        legend.box.background = element_rect(colour = "black"),
        legend.background = element_blank(),
        legend.key = element_blank(), 
        panel.background = element_rect(fill = 'white', colour = 'white'),
        panel.grid.major = element_line(colour = "grey95"),
        panel.grid.major.x = element_blank(), 
        panel.grid.minor = element_line(colour = "grey95"),
        panel.grid.minor.x = element_blank(), 
        axis.line.x = element_line(color="black"),
        axis.line.y = element_line(color="black"),
        strip.background = element_rect(fill="white")) 
fuji_plot
ggsave('03_figures/fg_02.jpg', 
       width = 9, 
       height = 6)

# Table A4 (Supplemental Material): Policy Support Regression Models -----------

# Linear Models (Coefficients)
m0_cf_c <- lm(position  ~ cond*cha_id, 
              data = df_e[df_e$exp_type == "fuji", ]) # Baseline
m0_ca_c <- lm(position  ~ cond*cha_id, 
              data = df_e[df_e$exp_type == "anti-fuji", ]) # Baseline
m1_cf_c <- lm(position  ~ cond*cha_id + party_id + ido_id + pb06 + order + exp, 
              data = df_e[df_e$exp_type == "fuji", ]) # controls
m1_ca_c <- lm(position ~ cond*cha_id + party_id + ido_id + pb06 + order + exp, 
              data = df_e[df_e$exp_type == "anti-fuji", ]) # controls
# Standard Errors
m0_cf <- lm.cluster(position  ~ cond*cha_id, 
                    data = df_e[df_e$exp_type == "fuji", ],
                    cluster = 'ResponseId') # Baseline
m0_ca <- lm.cluster(position  ~ cond*cha_id, 
                    data = df_e[df_e$exp_type == "anti-fuji", ],
                    cluster = 'ResponseId') # Baseline
m1_cf <- lm.cluster(position  ~ cond*cha_id + party_id + ido_id + pb06 + order + exp, 
                    data = df_e[df_e$exp_type == "fuji", ],
                    cluster = 'ResponseId') # controls
m1_ca <- lm.cluster(position  ~ cond*cha_id + party_id + ido_id + pb06 + order + exp, 
                    data = df_e[df_e$exp_type == "anti-fuji", ],
                    cluster = 'ResponseId') # controls
# Hierarchical Models
m0_mf <- lmer(position ~ cond*cha_id +
                (1 | ResponseId), data = df_e[df_e$exp_type == "fuji", ]) # Baseline
m0_ma <- lmer(position ~ cond*cha_id + 
                (1 | ResponseId), data = df_e[df_e$exp_type == "anti-fuji", ]) # Baseline
m1_mf <- lmer(position ~ cond*cha_id + party_id + ido_id + pb06 + order + exp  + 
                (1 | ResponseId), data = df_e[df_e$exp_type == "fuji", ]) # Controls
m1_ma <- lmer(position ~ cond*cha_id + party_id + ido_id + pb06 + order + exp + 
                (1 | ResponseId), data = df_e[df_e$exp_type == "anti-fuji", ]) # Controls

# Table with main coefficients
stargazer(m0_cf_c,
          m0_mf,
          m1_cf_c,
          m1_mf,
          m0_ca_c,
          m0_ma,
          m1_ca_c,
          m1_ma,
          se = list(summary(m0_cf)[, 2],
                    coef(summary(m0_mf))[ , "Std. Error"],
                    summary(m1_cf)[, 2],
                    coef(summary(m1_mf))[ , "Std. Error"],
                    summary(m0_ca)[, 2],
                    coef(summary(m0_ma))[ , "Std. Error"],
                    summary(m1_ca)[, 2],
                    coef(summary(m1_ma))[ , "Std. Error"]), 
          # type = "text",
          type = "html", 
          out = "04_outputs/table_a04.doc",
          dep.var.labels = c('Policy Support'), 
          column.labels   = c("Fujimorista Model", "Anti-Fujimorista Model"),
          column.separate = c(4, 4),
          covariate.labels = c("Fujimorista Cue (ref. control)",
                               "Anti-Fujimorista Cue (ref. control)",
                               "Fujimorista (ref. none)",
                               "Anti-Fujimorista (ref. none)",
                               "Fujimorista Cue $\\cdot$ Fujimorista",
                               "Fujimorista Cue $\\cdot$ Anti-Fujimorista",
                               "Anti-Fujimorista Cue $\\cdot$ Fujimorista",
                               "Anti-Fujimorista Cue $\\cdot$ Anti-Fujimorista",
                               "Intercept"),
          omit = c("party_id", "ido_id", "pb06", "order", "exp"),
          omit.stat = c("ll", "aic", 'bic', 'f', 'ser'),
          add.lines = list(c("Conditional R$^{2}$", 
                             "",
                             round(as.numeric(performance::r2(m0_mf)[1]), 3),
                             "",
                             round(as.numeric(performance::r2(m1_mf)[1]), 3), 
                             "",
                             round(as.numeric(performance::r2(m0_ma)[1]), 3),
                             "",
                             round(as.numeric(performance::r2(m1_ma)[1]), 3)),
                           c("Marginal R$^{2}$", 
                             "",
                             round(as.numeric(performance::r2(m0_mf)[2]), 3),
                             "",
                             round(as.numeric(performance::r2(m1_mf)[2]), 3), 
                             "",
                             round(as.numeric(performance::r2(m0_ma)[2]), 3),
                             "",
                             round(as.numeric(performance::r2(m1_ma)[2]), 3)),
                           c("Model", "OLS", "HLM", "OLS", "HLM", "OLS", "HLM", "OLS", "HLM"),
                           c("Controls", "No", "No", "Yes", "Yes", "No", "No", "Yes", "Yes")),
          table.layout = "=dc-#-t-sa-n",
          star.cutoffs = c(0.05),
          notes        = "OLS models estimated with robust standard error at the respondent level."
          )

# Table with all coefficients
stargazer(m0_cf_c,
          m0_mf,
          m1_cf_c,
          m1_mf,
          m0_ca_c,
          m0_ma,
          m1_ca_c,
          m1_ma,
          se = list(summary(m0_cf)[, 2],
                    coef(summary(m0_mf))[ , "Std. Error"],
                    summary(m1_cf)[, 2],
                    coef(summary(m1_mf))[ , "Std. Error"],
                    summary(m0_ca)[, 2],
                    coef(summary(m0_ma))[ , "Std. Error"],
                    summary(m1_ca)[, 2],
                    coef(summary(m1_ma))[ , "Std. Error"]), 
          # type = "text",
          type = "html", 
          out = "04_outputs/table_a04_complete.doc",
          dep.var.labels = c('Policy Support'), 
          column.labels   = c("Fujimorista Model", "Anti-Fujimorista Model"),
          column.separate = c(4, 4),
          order = c(1:4, 20:21, 24:25, 5, 8, 7, 6, 9:11, 14, 15, 13, 16, 12),
          covariate.labels = c("Fujimorista Cue (ref. control)",
                               "Anti-Fujimorista Cue (ref. control)",
                               "Fujimorista (ref. none)",
                               "Anti-Fujimorista (ref. none)",
                               "Fujimorista Cue $\\cdot$ Fujimorista",
                               "Fujimorista Cue $\\cdot$ Anti-Fujimorista",
                               "Anti-Fujimorista Cue $\\cdot$ Fujimorista",
                               "Anti-Fujimorista Cue $\\cdot$ Anti-Fujimorista",
                               "Party ID: Fuerza Popular (ref. Acción Popular)", 
                               "Party ID: Perú Libre (ref. Acción Popular)",
                               "Party ID: Other (ref. Acción Popular)",
                               "Party ID: None (ref. Acción Popular)",
                               "Ideological ID: Center (ref. Left)",
                               "Ideological ID: Right (ref. Left)",
                               "Ideological ID: None (ref. Left)",
                               "Pol. Interest: Don't know (ref. Interested)",
                               "Pol. Interest: Not at all interested (ref. Interested)",
                               "Pol. Interest: Not very interested (ref. Interested)",
                               "Pol. Interest: Somewhat interested (ref. Interested)",
                               "Pol. Interest: Very interested (ref. Interested)",
                               "Experiment Order", 
                               "Experiment 4 (ref. experiment 1)", 
                               "Experiment 6 (ref. experiment 1)", 
                               "Experiment 3 (ref. experiment 2)", 
                               "Experiment 5 (ref. experiment 2)", 
                               "Intercept"),
          omit.stat = c("ll", "aic", 'bic', 'f', 'ser'),
          add.lines = list(c("Conditional R$^{2}$", 
                             "",
                             round(as.numeric(performance::r2(m0_mf)[1]), 3),
                             "",
                             round(as.numeric(performance::r2(m1_mf)[1]), 3), 
                             "",
                             round(as.numeric(performance::r2(m0_ma)[1]), 3),
                             "",
                             round(as.numeric(performance::r2(m1_ma)[1]), 3)),
                           c("Marginal R$^{2}$", 
                             "",
                             round(as.numeric(performance::r2(m0_mf)[2]), 3),
                             "",
                             round(as.numeric(performance::r2(m1_mf)[2]), 3), 
                             "",
                             round(as.numeric(performance::r2(m0_ma)[2]), 3),
                             "",
                             round(as.numeric(performance::r2(m1_ma)[2]), 3)),
                           c("Model", "OLS", "HLM", "OLS", "HLM", "OLS", "HLM", "OLS", "HLM"),
                           c("Controls", "No", "No", "Yes", "Yes", "No", "No", "Yes", "Yes")),
          table.layout = "=dc-#-t-sa-n",
          star.cutoffs = c(0.05),
          notes        = "OLS models estimated with robust standard error at the respondent level.")

# Fujimorista Cues plot --------------------------------------------------------

# Set seed and number of iterations
set.seed(20220411)

# FUNCTION
mm <- list(
  no = data.frame(rbind(
    c(0, 1, 0, 0, 1, 0),
    c(0, 1, 0, 0, 0, 1),
    c(0, 1, 0, 0, 0, 0)
  )),
  yes = data.frame(rbind(
    c(0, 1, 0, 0, rep(0, 15), 1, 0),
    c(0, 1, 0, 0, rep(0, 15), 0, 1),
    c(0, 1, 0, 0, rep(0, 15), 0, 0)
  ))
)

sim_preds <- function(iterations, model, model_type = "lm", controls = "no") {
  # Beta Simulations
  beta.sim <- MASS::mvrnorm(iterations, 
                              mu = if(model_type == "lm") coef(model) else coef(summary(model))[,1], 
                              Sigma = vcov(model)) # variance of the Betas
  # Predictions
  preds <- as.matrix(mm[[controls]]) %*% t(as.matrix(beta.sim))
  # Quantiles
  preds_q <- t(apply(preds, 1, quantile, probs = c(0.025, 0.5, 0.975)))
  # Variables
  preds_v <- cbind(preds_q,
                  mm[[controls]][, (ncol(mm[[controls]]) - 1):ncol(mm[[controls]])])
  names(preds_v) <- c("lwr", "scr", "upr", "fuji", "anti")
  preds_v$model <- model_type
  preds_v$controls <- controls
  return(preds_v)
}

# Predictions
preds <- bind_rows(sim_preds(5000, m0_cf, "lm", "no"), 
                   sim_preds(5000, m1_cf, "lm", "yes"),
                   sim_preds(5000, m0_mf, "hlm", "no"),
                   sim_preds(5000, m1_mf, "hlm", "yes"))

# PLOT
fuji_cue <- preds %>% 
  mutate(cha = ifelse(fuji == 1, "Fujimorista", 
                      ifelse(anti == 1, "Anti-Fujimorista", 
                             "None")), 
         cha = factor(cha, levels = c("Fujimorista", 
                                      "None", 
                                      "Anti-Fujimorista")), 
         mo_c = case_when(
           model == "lm" & controls == "no" ~ "LM",
           model == "lm" & controls == "yes" ~ "LM-Controls",
           model == "hlm" & controls == "no" ~ "HLM",
           model == "hlm" & controls == "yes" ~ "HLM-Controls"), 
         mo_c = factor(mo_c, levels = c("LM", 
                                        "HLM", 
                                        "LM-Controls", 
                                        "HLM-Controls"))
  ) %>% 
  ggplot(aes(x = cha, y = scr, colour = mo_c, shape = mo_c)) +
  geom_point(position = position_dodge(width = 0.5)) +
  geom_errorbar(aes(ymin = lwr, ymax = upr), 
                width = 0, 
                position = position_dodge(width = 0.5)) +
  scale_shape_manual(values = c(16, 16, 17, 17)) +
  scale_color_manual(values = c("black", "grey", "black", "grey")) +
  scale_y_continuous(limits = c(-1.5, 1.5)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = 'Identity', y = 'Marginal Effect', title = "Fujimorista Cue") +
  theme(legend.position = "bottom", 
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5),
        legend.box.background = element_rect(colour = "black"),
        legend.background = element_blank(),
        legend.key = element_blank(),
        panel.background = element_rect(fill = 'white', colour = 'white'),
        panel.grid.major = element_line(colour = "grey95"),
        panel.grid.major.x = element_blank(),
        panel.grid.minor = element_line(colour = "grey95"),
        panel.grid.minor.x = element_blank(),
        axis.line.x = element_line(color="black"),
        axis.line.y = element_line(color="black"),
        strip.background = element_rect(fill="white")) +
  guides(colour = guide_legend(nrow = 2), 
         shape = guide_legend(nrow = 2))
fuji_cue

# Anti- Fujimorista Cues plot --------------------------------------------------

# Predictions
preds <- bind_rows(sim_preds(5000, m0_ca, "lm", "no"), 
                   sim_preds(5000, m1_ca, "lm", "yes"),
                   sim_preds(5000, m0_ma, "hlm", "no"),
                   sim_preds(5000, m1_ma, "hlm", "yes"))

# PLOT
anti_cue <- preds %>% 
  mutate(cha = ifelse(fuji == 1, "Fujimorista", 
                      ifelse(anti == 1, "Anti-Fujimorista", 
                             "None")), 
         cha = factor(cha, levels = c("Fujimorista", 
                                      "None", 
                                      "Anti-Fujimorista")), 
         mo_c = case_when(
           model == "lm" & controls == "no" ~ "LM",
           model == "lm" & controls == "yes" ~ "LM-Controls",
           model == "hlm" & controls == "no" ~ "HLM",
           model == "hlm" & controls == "yes" ~ "HLM-Controls"), 
         mo_c = factor(mo_c, levels = c("LM", 
                                        "HLM", 
                                        "LM-Controls", 
                                        "HLM-Controls"))
  ) %>% 
  ggplot(aes(x = cha, y = scr, colour = mo_c, shape = mo_c)) +
  geom_point(position = position_dodge(width = 0.5)) +
  geom_errorbar(aes(ymin = lwr, ymax = upr), 
                width = 0, 
                position = position_dodge(width = 0.5)) +
  scale_shape_manual(values = c(16, 16, 17, 17)) +
  scale_color_manual(values = c("black", "grey", "black", "grey")) +
  scale_y_continuous(limits = c(-1.5, 1.5)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = 'Identity', y = 'Marginal Effect', title = "Anti-Fujimorista Cue") +
  theme(legend.position = "bottom", 
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5),
        legend.box.background = element_rect(colour = "black"),
        legend.background = element_blank(),
        legend.key = element_blank(), 
        panel.background = element_rect(fill = 'white', colour = 'white'),
        panel.grid.major = element_line(colour = "grey95"),
        panel.grid.major.x = element_blank(),
        panel.grid.minor = element_line(colour = "grey95"),
        panel.grid.minor.x = element_blank(), 
        axis.line.x = element_line(color="black"),
        axis.line.y = element_line(color="black"),
        strip.background = element_rect(fill="white")) +
  guides(colour = guide_legend(nrow = 2), 
         shape = guide_legend(nrow = 2))
anti_cue

# Figure 4: Fujimorista and Anti-Fujimorista Cues on identifiers (discrete) ----

# Merge plots
ggarrange(fuji_cue, 
          anti_cue, 
          common.legend = TRUE, 
          legend = 'bottom') +
  theme(plot.margin = margin(b = 1))

ggsave("03_figures/fg_04.jpg", 
       width = 12, 
       height = 6)

# Figure A1: Individual Cues ---------------------------------------------------

# BASE MODELS
m1_1e <- lm(position  ~ cond*cha_id, 
            data = df_e[df_e$exp == 1, ])
m1_2e <- lm(position  ~ cond*cha_id, 
            data = df_e[df_e$exp == 2, ])
m1_3e <- lm(position  ~ cond*cha_id, 
            data = df_e[df_e$exp == 3, ])
m1_4e <- lm(position  ~ cond*cha_id, 
            data = df_e[df_e$exp == 4, ])
m1_5e <- lm(position  ~ cond*cha_id, 
            data = df_e[df_e$exp == 5, ])
m1_6e <- lm(position  ~ cond*cha_id, 
            data = df_e[df_e$exp == 6, ])

pred_1e <- sim_preds(5000, m1_1e) %>% mutate(exp = "Experiment 1")
pred_2e <- sim_preds(5000, m1_2e) %>% mutate(exp = "Experiment 2")
pred_3e <- sim_preds(5000, m1_3e) %>% mutate(exp = "Experiment 3")
pred_4e <- sim_preds(5000, m1_4e) %>% mutate(exp = "Experiment 4")
pred_5e <- sim_preds(5000, m1_5e) %>% mutate(exp = "Experiment 5")
pred_6e <- sim_preds(5000, m1_6e) %>% mutate(exp = "Experiment 6")

preds_e <- rbind(pred_1e, 
                 pred_2e,
                 pred_3e,
                 pred_4e,
                 pred_5e,
                 pred_6e) %>% 
  mutate(cha = ifelse(fuji == 1, "Fujimorista", 
                      ifelse(anti == 1, "Anti-Fujimorista", 
                             "None")), 
         cha = factor(cha, levels = c("Fujimorista", 
                                      "None", 
                                      "Anti-Fujimorista")))

preds_e %>% ggplot(aes(x = cha, y = scr)) +
  geom_point(position = position_dodge(width = 0.5)) +
  geom_errorbar(aes(ymin = lwr, ymax = upr), 
                width = 0, 
                position = position_dodge(width = 0.5)) +
  scale_y_continuous(limits = c(-2.5, 2.5)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = 'Identity', y = 'Marginal Effect') +
  facet_wrap(~ exp) +
  theme(strip.background = element_rect(fill = "white", colour = "black"),
        panel.background = element_rect(fill = 'white', colour = 'black'),
        panel.grid.major = element_line(colour = "grey95"),
        panel.grid.major.x = element_blank(), # Remove vertical lines
        panel.grid.minor = element_line(colour = "grey95"),
        panel.grid.minor.x = element_blank())
ggsave('03_figures/fg_A01.jpg', 
       width = 10, 
       height = 6)

# Table A5 (Supplemental Material): Pol. Support Regressions M. (strength) -----

# Linear Models (Coefficients)
m0_cf_c <- lm(position  ~ cond*fuji_st + cond*anti_st, 
              data = df_e[df_e$exp_type == "fuji", ]) # Baseline
m0_ca_c <- lm(position  ~ cond*fuji_st + cond*anti_st, 
              data = df_e[df_e$exp_type == "anti-fuji", ]) # Baseline
m1_cf_c <- lm(position  ~ cond*fuji_st + cond*anti_st + party_id + ido_id + pb06 + order + exp, 
              data = df_e[df_e$exp_type == "fuji", ]) # controls
m1_ca_c <- lm(position ~ cond*fuji_st + cond*anti_st + party_id + ido_id + pb06 + order + exp, 
              data = df_e[df_e$exp_type == "anti-fuji", ]) # controls
# Standard Errors
m0_cf <- lm.cluster(position  ~ cond*fuji_st + cond*anti_st, 
                    data = df_e[df_e$exp_type == "fuji", ],
                    cluster = 'ResponseId') # Baseline
m0_ca <- lm.cluster(position  ~ cond*fuji_st + cond*anti_st, 
                    data = df_e[df_e$exp_type == "anti-fuji", ],
                    cluster = 'ResponseId') # Baseline
m1_cf <- lm.cluster(position  ~ cond*fuji_st + cond*anti_st + party_id + ido_id + pb06 + order + exp, 
                    data = df_e[df_e$exp_type == "fuji", ],
                    cluster = 'ResponseId') # controls
m1_ca <- lm.cluster(position  ~ cond*fuji_st + cond*anti_st + party_id + ido_id + pb06 + order + exp, 
                    data = df_e[df_e$exp_type == "anti-fuji", ],
                    cluster = 'ResponseId') # controls
# Hierarchical Models
m0_mf <- lmer(position ~ cond*fuji_st + cond*anti_st +
                (1 | ResponseId), data = df_e[df_e$exp_type == "fuji", ]) # Baseline
m0_ma <- lmer(position ~ cond*fuji_st + cond*anti_st + 
                (1 | ResponseId), data = df_e[df_e$exp_type == "anti-fuji", ]) # Baseline
m1_mf <- lmer(position ~ cond*fuji_st + cond*anti_st + party_id + ido_id + pb06 + order + exp  + 
                (1 | ResponseId), data = df_e[df_e$exp_type == "fuji", ]) # Controls
m1_ma <- lmer(position ~ cond*fuji_st + cond*anti_st + party_id + ido_id + pb06 + order + exp + 
                (1 | ResponseId), data = df_e[df_e$exp_type == "anti-fuji", ]) # Controls
# Table with Main coefficients
stargazer(m0_cf_c,
          m0_mf,
          m1_cf_c,
          m1_mf,
          m0_ca_c,
          m0_ma,
          m1_ca_c,
          m1_ma,
          se = list(summary(m0_cf)[, 2],
                    coef(summary(m0_mf))[ , "Std. Error"],
                    summary(m1_cf)[, 2],
                    coef(summary(m1_mf))[ , "Std. Error"],
                    summary(m0_ca)[, 2],
                    coef(summary(m0_ma))[ , "Std. Error"],
                    summary(m1_ca)[, 2],
                    coef(summary(m1_ma))[ , "Std. Error"]), 
          # type = "text",
          type = "html", 
          out = "04_outputs/table_a05.doc",
          dep.var.labels = c('Policy Support'), 
          column.labels   = c("Fujimorista Model", "Anti-Fujimorista Model"),
          column.separate = c(4, 4),
          covariate.labels = c("Fujimorista Cue (ref. control)",
                               "Anti-Fujimorista Cue (ref. control)",
                               "Fujimorista Identity Strength",
                               "Anti-Fujimorista Identity Strength",
                               "Fujimorista Cue $\\cdot$ Fujimorista",
                               "Fujimorista Cue $\\cdot$ Anti-Fujimorista",
                               "Anti-Fujimorista Cue $\\cdot$ Fujimorista",
                               "Anti-Fujimorista Cue $\\cdot$ Anti-Fujimorista",
                               "Intercept"),
          omit = c("party_id", "ido_id", "pb06", "order", "exp"),
          omit.stat = c("ll", "aic", 'bic', 'f', 'ser'),
          add.lines = list(c("Conditional R$^{2}$", 
                             "",
                             round(as.numeric(performance::r2(m0_mf)[1]), 3),
                             "",
                             round(as.numeric(performance::r2(m1_mf)[1]), 3), 
                             "",
                             round(as.numeric(performance::r2(m0_ma)[1]), 3),
                             "",
                             round(as.numeric(performance::r2(m1_ma)[1]), 3)),
                           c("Marginal R$^{2}$", 
                             "",
                             round(as.numeric(performance::r2(m0_mf)[2]), 3),
                             "",
                             round(as.numeric(performance::r2(m1_mf)[2]), 3), 
                             "",
                             round(as.numeric(performance::r2(m0_ma)[2]), 3),
                             "",
                             round(as.numeric(performance::r2(m1_ma)[2]), 3)),
                           c("Model", "OLS", "HLM", "OLS", "HLM", "OLS", "HLM", "OLS", "HLM"),
                           c("Controls", "No", "No", "Yes", "Yes", "No", "No", "Yes", "Yes")),
          table.layout = "=dc-#-t-sa-n",
          star.cutoffs = c(0.05),
          notes        = "OLS models estimated with robust standard error at the respondent level."
)

# Table complete
stargazer(m0_cf_c,
          m0_mf,
          m1_cf_c,
          m1_mf,
          m0_ca_c,
          m0_ma,
          m1_ca_c,
          m1_ma,
          se = list(summary(m0_cf)[, 2],
                    coef(summary(m0_mf))[ , "Std. Error"],
                    summary(m1_cf)[, 2],
                    coef(summary(m1_mf))[ , "Std. Error"],
                    summary(m0_ca)[, 2],
                    coef(summary(m0_ma))[ , "Std. Error"],
                    summary(m1_ca)[, 2],
                    coef(summary(m1_ma))[ , "Std. Error"]), 
          # type = "text",
          type = "html",
          out = "04_outputs/table_a05_complete.doc",
          dep.var.labels = c('Policy Support'), 
          column.labels   = c("Fujimorista Model", "Anti-Fujimorista Model"),
          column.separate = c(4, 4),
          order = c(1:4, 20:21, 24:25, 5, 8, 7, 6, 9:11, 14, 15, 13, 16, 12),
          covariate.labels = c("Fujimorista Cue (ref. control)",
                               "Anti-Fujimorista Cue (ref. control)",
                               "Fujimorista Identity Strength",
                               "Anti-Fujimorista Identity Strength",
                               "Fujimorista Cue $\\cdot$ Fujimorista",
                               "Fujimorista Cue $\\cdot$ Anti-Fujimorista",
                               "Anti-Fujimorista Cue $\\cdot$ Fujimorista",
                               "Anti-Fujimorista Cue $\\cdot$ Anti-Fujimorista",
                               "Party ID: Fuerza Popular (ref. Acción Popular)",
                               "Party ID: Perú Libre (ref. Acción Popular)",
                               "Party ID: Other (ref. Acción Popular)",
                               "Party ID: None (ref. Acción Popular)",
                               "Ideological ID: Center (ref. Left)",
                               "Ideological ID: Right (ref. Left)",
                               "Ideological ID: None (ref. Left)",
                               "Pol. Interest: Don't know (ref. Interested)",
                               "Pol. Interest: Not at all interested (ref. Interested)",
                               "Pol. Interest: Not very interested (ref. Interested)",
                               "Pol. Interest: Somewhat interested (ref. Interested)",
                               "Pol. Interest: Very interested (ref. Interested)",
                               "Experiment Order",
                               "Experiment 4 (ref. experiment 1)",
                               "Experiment 6 (ref. experiment 1)",
                               "Experiment 3 (ref. experiment 2)",
                               "Experiment 5 (ref. experiment 2)",
                               "Intercept"),
          omit.stat = c("ll", "aic", 'bic', 'f', 'ser'),
          add.lines = list(c("Conditional R$^{2}$", 
                             "",
                             round(as.numeric(performance::r2(m0_mf)[1]), 3),
                             "",
                             round(as.numeric(performance::r2(m1_mf)[1]), 3), 
                             "",
                             round(as.numeric(performance::r2(m0_ma)[1]), 3),
                             "",
                             round(as.numeric(performance::r2(m1_ma)[1]), 3)),
                           c("Marginal R$^{2}$", 
                             "",
                             round(as.numeric(performance::r2(m0_mf)[2]), 3),
                             "",
                             round(as.numeric(performance::r2(m1_mf)[2]), 3), 
                             "",
                             round(as.numeric(performance::r2(m0_ma)[2]), 3),
                             "",
                             round(as.numeric(performance::r2(m1_ma)[2]), 3)),
                           c("Model", "OLS", "HLM", "OLS", "HLM", "OLS", "HLM", "OLS", "HLM"),
                           c("Controls", "No", "No", "Yes", "Yes", "No", "No", "Yes", "Yes")),
          table.layout = "=dc-#-t-sa-n",
          star.cutoffs = c(0.05),
          notes        = "OLS models estimated with robust standard error at the respondent level."
)

# Figure 5: Fujimorista and Anti-Fujimorista Cues (strength) -------------------

# Set seed and number of iterations
set.seed(20220411)

# FUNCTION
mm <- list(
  no = data.frame(rbind(
    c(0, 1, 0, 0, 0.0, 0),
    c(0, 1, 0, 0, 0.1, 0),
    c(0, 1, 0, 0, 0.2, 0),
    c(0, 1, 0, 0, 0.3, 0),
    c(0, 1, 0, 0, 0.4, 0),
    c(0, 1, 0, 0, 0.5, 0),
    c(0, 1, 0, 0, 0.6, 0),
    c(0, 1, 0, 0, 0.7, 0),
    c(0, 1, 0, 0, 0.8, 0),
    c(0, 1, 0, 0, 0.9, 0),
    c(0, 1, 0, 0, 1.0, 0),
    c(0, 1, 0, 0, 0, 0.0),
    c(0, 1, 0, 0, 0, 0.1),
    c(0, 1, 0, 0, 0, 0.2),
    c(0, 1, 0, 0, 0, 0.3),
    c(0, 1, 0, 0, 0, 0.4),
    c(0, 1, 0, 0, 0, 0.5),
    c(0, 1, 0, 0, 0, 0.6),
    c(0, 1, 0, 0, 0, 0.7),
    c(0, 1, 0, 0, 0, 0.8),
    c(0, 1, 0, 0, 0, 0.9),
    c(0, 1, 0, 0, 0, 1.0)
  )),
  yes = data.frame(rbind(
    c(0, 1, 0, 0, rep(0, 15), 0.001, 0),
    c(0, 1, 0, 0, rep(0, 15), 0.1, 0),
    c(0, 1, 0, 0, rep(0, 15), 0.2, 0),
    c(0, 1, 0, 0, rep(0, 15), 0.3, 0),
    c(0, 1, 0, 0, rep(0, 15), 0.4, 0),
    c(0, 1, 0, 0, rep(0, 15), 0.5, 0),
    c(0, 1, 0, 0, rep(0, 15), 0.6, 0),
    c(0, 1, 0, 0, rep(0, 15), 0.7, 0),
    c(0, 1, 0, 0, rep(0, 15), 0.8, 0),
    c(0, 1, 0, 0, rep(0, 15), 0.9, 0),
    c(0, 1, 0, 0, rep(0, 15), 1.0, 0),
    c(0, 1, 0, 0, rep(0, 15), 0, 0.001),
    c(0, 1, 0, 0, rep(0, 15), 0, 0.1),
    c(0, 1, 0, 0, rep(0, 15), 0, 0.2),
    c(0, 1, 0, 0, rep(0, 15), 0, 0.3),
    c(0, 1, 0, 0, rep(0, 15), 0, 0.4),
    c(0, 1, 0, 0, rep(0, 15), 0, 0.5),
    c(0, 1, 0, 0, rep(0, 15), 0, 0.6),
    c(0, 1, 0, 0, rep(0, 15), 0, 0.7),
    c(0, 1, 0, 0, rep(0, 15), 0, 0.8),
    c(0, 1, 0, 0, rep(0, 15), 0, 0.9),
    c(0, 1, 0, 0, rep(0, 15), 0, 1.0)
  ))
)

# Predictions
preds <- bind_rows(sim_preds(5000, m0_cf, "lm", "no"), 
                   sim_preds(5000, m1_cf, "lm", "yes"),
                   sim_preds(5000, m0_mf, "hlm", "no"),
                   sim_preds(5000, m1_mf, "hlm", "yes"))
# FUJI
fuji_st <- preds %>% 
  mutate(mo_c = case_when(
    model == "lm" & controls == "no" ~ "LM",
    model == "lm" & controls == "yes" ~ "LM-Controls",
    model == "hlm" & controls == "no" ~ "HLM",
    model == "hlm" & controls == "yes" ~ "HLM-Controls"), 
    mo_c = factor(mo_c, levels = c("LM", 
                                   "HLM", 
                                   "LM-Controls", 
                                   "HLM-Controls"))) %>% 
  pivot_longer(cols = c(fuji, anti)) %>% 
  mutate(name = ifelse(name == "fuji", "Fujimorista", "Anti-Fujimorista"), 
         name = factor(name, levels = c("Fujimorista", 
                                          "Anti-Fujimorista"))) %>% 
  filter(value != 0) %>%
  ggplot(aes(x = value, y = scr, colour = name, shape = mo_c)) +
  geom_point(position = position_dodge(width = 0.08)) +
  geom_errorbar(aes(ymin = lwr, ymax = upr), 
                width = 0, 
                position = position_dodge(width = 0.08)) +
  scale_shape_manual(values = c(15, 16, 17, 18)) +
  scale_color_manual(values = c("black", "grey", "black", "grey")) +
  scale_y_continuous(limits = c(-2.2, 2.2)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = 'Identity Strength', y = 'Marginal Effect', title = "Fujimorista Cue") +
  theme(legend.position = "bottom", 
        plot.title = element_text(hjust = 0.5),
        legend.box.background = element_rect(colour = "black"),
        legend.background = element_blank(),
        legend.key = element_blank(), 
        panel.background = element_rect(fill = 'white', colour = 'white'),
        panel.grid.major = element_line(colour = "grey95"),
        panel.grid.major.x = element_blank(), 
        panel.grid.minor = element_line(colour = "grey95"),
        panel.grid.minor.x = element_blank(), 
        axis.line.x = element_line(color="black"),
        axis.line.y = element_line(color="black"),
        strip.background = element_rect(fill="white")) +
  guides(colour = guide_legend("Identity",
                               title.position = "top", 
                               title.hjust = 0.5,
                               nrow = 2),
         shape = guide_legend("Model",
                              title.position = "top", 
                              title.hjust = 0.5, 
                              nrow = 2))
fuji_st

# Predictions
preds <- bind_rows(sim_preds(5000, m0_ca, "lm", "no"), 
                   sim_preds(5000, m1_ca, "lm", "yes"),
                   sim_preds(5000, m0_ma, "hlm", "no"),
                   sim_preds(5000, m1_ma, "hlm", "yes"))
# ANTI
anti_st <- preds %>% 
  mutate(mo_c = case_when(
    model == "lm" & controls == "no" ~ "LM",
    model == "lm" & controls == "yes" ~ "LM-Controls",
    model == "hlm" & controls == "no" ~ "HLM",
    model == "hlm" & controls == "yes" ~ "HLM-Controls"), 
    mo_c = factor(mo_c, levels = c("LM", 
                                   "HLM", 
                                   "LM-Controls", 
                                   "HLM-Controls"))) %>% 
  pivot_longer(cols = c(fuji, anti)) %>% 
  mutate(name = ifelse(name == "fuji", "Fujimorista", "Anti-Fujimorista"), 
         name = factor(name, levels = c("Fujimorista", 
                                        "Anti-Fujimorista"))) %>% 
  filter(value != 0) %>% 
  ggplot(aes(x = value, y = scr, colour = name, shape = mo_c)) +
  geom_point(position = position_dodge(width = 0.08)) +
  geom_errorbar(aes(ymin = lwr, ymax = upr), 
                width = 0, 
                position = position_dodge(width = 0.08)) +
  scale_shape_manual(values = c(15, 16, 17, 18)) +
  scale_color_manual(values = c("black", "grey", "black", "grey")) +
  scale_y_continuous(limits = c(-2.2, 2.2)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = 'Identity Strength', y = 'Marginal Effect', title = "Anti-Fujimorista Cue") +
  theme(legend.position = "bottom", 
        plot.title = element_text(hjust = 0.5),
        legend.box.background = element_rect(colour = "black"),
        legend.background = element_blank(), 
        legend.key = element_blank(),
        panel.background = element_rect(fill = 'white', colour = 'white'),
        panel.grid.major = element_line(colour = "grey95"),
        panel.grid.major.x = element_blank(), 
        panel.grid.minor = element_line(colour = "grey95"),
        panel.grid.minor.x = element_blank(),
        axis.line.x = element_line(color="black"),
        axis.line.y = element_line(color="black"),
        strip.background = element_rect(fill="white")) +
  guides(colour = guide_legend("Identity",
                               title.position = "top", 
                               title.hjust = 0.5,
                               nrow = 2),
         shape = guide_legend("Model",
                              title.position = "top", 
                              title.hjust = 0.5, 
                              nrow = 2))
anti_st

# All together
ggarrange(fuji_st, 
          anti_st, 
          common.legend = TRUE, 
          legend = "bottom") +
  theme(plot.margin = margin(b = 1))

ggsave("03_figures/fg_05.jpg", 
       width = 12, 
       height = 6)